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CHAPTER  I  • 

THE  PROBLEM 

INTRODUCTION 

This  paper  is  concerned  with  the  following  type  of  optimal  control  prob- 

lem: 

For  a  given  topological  configuration  of  a  linear  control 
system  and  a  given  input  f(t),  the  parameter  values  are  to 
be  found  which  optimize  some  performance  measure. 

Traditionally  problems  of  this  type  have  been  solved  by  trial  and  error 
(Bach,  1965).   However  this  is  time  consuming,  costly,  and  impractical  on 
large  systems.   In  recent  years  considerable  interest  has  developed  in  using 
computers  as  a  design  tool  to  simulate  the  control  system  and  find  the  most 
desirable  system  parameters.   This  is  done  by  defining  a  performance  measure 
or  measures  in  terms  of  a  distance  function.   A  distance  function  relates 
quality  of  performance  to  the  least  distance  from  the  origin.   Any  number  of 
performance  measures  can  be  combined  in  a  single  distance  function.   An  exam- 
ple of  a  commonly  used  distance  function  is  the  square  root  of  the  sum  of  the 
squares  of  the  performance  criterion  (Levine,  1964). 

In  control  systems  the  performance  measure  usually  involves  the  integral 
of  some  form  of  error  between  the  input  and  the  output.   Other  criterion  that 
may  be  used  are  overshoot,  undershoot,  time  delay,  and  settling  time  (Kuo, 
1962).   These  measures  are  defined  in  Appendix  A.  " 

Most  other  treatments  of  this  problem  are  restricted  to  step  inputs. 
This  approach  allows  a  more  general  f(t)  that  has  a  constant  final  value  after 


some  time  T^ .  The  assumption  is  made  in  this  paper  that  the  input  f (t)  never 
exceeds  the  final  value.   If  the  input  f(t)  is  larger  than  the  final  value  for 
some  time  T_  less  than  T^ ,  then  the  performance  measures  must  be  redefined  to 
'    take  this  into  account.  ' 

This  paper  describes  a  general  FORTRAN  program  to  solve  the  above  prob-  '- 
lem.   The  program  is  written  to  find  the  optimal  parameters  of  any  control 
system  that  can  be  described  by  an  W^   order  (2<N<8)  linear  differential  eq- 
uation with  constant  coefficients.   The  computation  time  increases  very  rap- 
idly however,  as  the  order  of  the  system  increases.   To  demonstrate  the  prog- 
ram, a  specific  problem  is  considered  and  the  performance  measure  criterion 
are  evaluated  for  two  distance  functions  and  three  different  error  definitions. 

STATEMENT  OF  THE  PROBLEM 

The  program  to  be  solved  is  as  follows.   The  parameters  of  a  third  order 
linear  differential  equation  are  to  be  found  which  minimize  a  distance  func- 
tion involving  settling  time  T  and  the  integral  of  some  form  of  error.   The 
input  to  the  system  is  a  modified  step  function  consisting  of  a  unit  ramp  to 
time  one  and  unit  step  from  one  to  infinity.   Thus  the  system  has  the  follow- 
ing form. 

•y  +  A^y   +  A^y   +   y   =   f(t) 
f(t)=0         -(»>t>Q 
t         0   ^  t  >  1 
1         1   ^  t 

The  distance  function  is  the  square  root  of  the  sum  of  the  squares  of  T 

s 

and  ERROR  where  T   is  defined  as  the  time  required  for  the  response  to  de- 
crease to  and  stay  within  five  percent  of  the  final  value.   Three  different 


ERROR  definitions  are  used  and  compared. 


1.    lAE  =    /  ERR  dt 


Integral  of  the  absolute  value  of  error. 

Note:  ERR  is  the  difference  between  the  input  and  output, 

T  can  either  be  a  fixed  value  designated  by  the 

user  or  the  value  of  T  . 

s 

T 

2.  lES  =    f   ERR^dt 

0 

Integral  of  the  error  squared. 

T 

3.  ITES=    /  t*ERR^dt 

0 

Integral  of  time  multiplied  error  squared. 


CHAPTER  II 
THEORETICAL  ANALYSIS, 

DISCUSSION 

Because  of  the  difficulties  in  generalizing  analytical  methods  of  opti- 
mization, search  techniques  are  used  to  find  the  optimal  parameter  values. 
The  distance  function  is  evaluated  each  time  the  parameters  are  changed  by 
solving  the  differential  equation  that  describes  the  system.   The  distance 
function  is  then  minimized  by  adjusting  the  parameters  in  some  optimal  way. 

For  efficiency  two  search  methods,  steepest  descent  and  relaxation,  are 
used.   Steepest  descent  works  best  in  adjusting  the  parameters  from  the  ini- 
tial values  to  somewhere  near  the  optimum  values.  This  method  has  problems 
with  step  sizes  near  the  minimum  and  doesn't  work  when  there  are  discontin- 
uities in  the  solution  space. 

Due  to  the  definition  of  T  ,  the  solution  space  for  the  distance  func- 

s' 

tion  is  discontinuous.   Near  the  optimal  parameter  values,  a  slight  change 
will  move  the  solution  outside  the  allowed  range  and  markedly  increase  the 
value  of  T  .   This  in  turn  increases  the  ERR  and  because  it  is  integrated 
over  a  longer  time. 

In  this  problem  the  minimum  values  occur  right  next  to  a  discontinuity 
in  the  solution  space.   (See  Fig.  1)   This  is  quite  understandable  when  one 
looks  at  the  definition  of  the  distance  function.   The  program  will  first  try 
and  make  the  response  converge  as  soon  as  possible  to  minimize  both  T  and 
ERR.   Then  once  T   is  relatively  constant,  the  program  will  move  the  overshoot 
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as  close  to  the  Limit  as  possible  without  exceeding  it.  Thus  it  can  be  seen 
why  a  very  small  increment  will  suddenly  cause  a  great  change  in  the  value  of 
the  distance  function.   This  causes  problems  in  the  search  technique  because 
proper  adjustments  must  be  made  when  this  happens.   It  is  also  easy  to  get 
caught  on  a  relative  minimum  near  a  discontinuity. 

This  brings  up  questions  as  to  the  validity  of  the  choice  of  performance 
measure.   One  of  the  requirements  for  an  optimum  system  is  mathematical  trac- 
tability  (Hancock,  1966).   Obviously  the  choice  of  the  performance  measure  is 
purely  subjective.   Therefore  to  be  effective  it  must  be  based  on  good  engi- 
neering judgement.   This  particular  performance  measure  should  only  be  used 
when  necessary.   In  most  problems,  other  measures  that  have  continuous  solu- 
tion spaces  would  be  more  satisfactory. 

As  noted  above  the  optimum  parameter  values  represent  a  relative  minimum 
and  not  necessarily  an  absolute  minimum.   When  there  are  several  minima  it  is 
necessary  to  run  the  program  several  times  from  different  starting  points  to 
determine  the  absolute  minimum. 

Because  the  parameters  represent  physical  components  the  accuracy  is  lim- 
ited to  one  part  in  a  thousand.   It  would  be  senseless  to  obtain  optimal  para- 
meter values  with  greater  accuracy  than  this  because  the  components  would 
change  with  time  or  temperature.   The  parameters  should  not  be  picked  too  near 
a  discontinuity  in  the  solution  space  for  the  same  reason. 

This  paper  assumes  there  should  be  no  steady  state  error.   This  con- 
strains the  parameter  on  y  to  the  value  1.   Thus  for  an  ISF''  order  differential 
equation  there  will  be  N-1  adjustable  parameters. 

Analysis  of  Methods  Used 

Three  numerical  methods  are  used  in  this  paper.   Two  are  search  tech- 


niquGS  and  the  other  is  a  method  for  finding  the  solution  to  a  differential 
equation.   The  methods  are  reviewed  briefly  here. 

Relaxation 

In  this  search  technique  a  single  parameter  is  varied  with  all  others 
held  constant  until  a  local  minimum  is  reached.   Slope  is  used  to  determine 
the  direction  of  descent.   Each  time  the  slope  changes,  the  stepsize  is  re- 
duced by  an  order  of  magnitude.   This  continues  until  some  termination  crit- 
erion is  satisfied.   When  the  solution  space  is  discontinuous  some  additional 
checks  are  made  to  insure  convergence. 

Steepest  Descent 

The  distance  function  is  defined  as  B(k^, . . . ,A^) .      This  function  can  be 
minimized  in  a  region  R  of  the  N-  dimensional  rectilinear  parameter  space  if 
one  starts  from  an  arbitrarily  selected  point  in  the  region  and  follows  a 
path  which  leads  to  decreasing  values  of  D.   The  rate  of.  change  dD/dt  is 
greatest  if  the  path  chosen  in  the  parameter  space  is  tangential  to  the  grad- 
ient vector  in  this  same  space  (Levine,  1964).   This  may  be  written  as  an  in- 
ner product  as  follows:         dD/dt  =  (grad  D)  •  (A) 

where  A  =  A^u^+A  u  +  ...  A  ii 

—    1—1   z— z       n— n 

The  u.'s  are  unit  vectors  along  the  coordinate  axis.   The  minimum  occurs  when 

the  two  vectors  are  parallel  and  point  in  opposite  directions  (Levine,  1964). 

This  may  be  written  as  follows: 

dA./dt  =  -(5D/^A.)    i  =  l,2,...n 

This  may  be  approximated  by  the  following  difference  equation: 

AA.  =  -(AD/AA,)At     i  =  1,2, ...n 

One  variation  on  steepest  descent  that  saves  some  computer  time  is  to 


calculate  the  steepest  descent  stepsizes  at  a  point  and  increment  with  those 
same  stepsizes  until  a  local  minimum  is  reached.  The  local  minimum  is  taken 
as  the  new  starting  point  and  the  process  is  repeated. 

Euler's  Method 

Euler's  Method  is  used  to  solve  the  differential  equation  because  of  its 
simplicity  and  minimal  computation  time.   If  higher  accuracy  would  be  needed 
some  other  method  such  as  Runge-Kutta  could  be  used  (Conte,  1965).   The 
choice  of  the  time  increment  has  an  important  part  in  the  accuracy  of  the  sol- 
ution, the  accuracy  of  the  settling  time  T  ,  and  the  computation  time  re- 
quired to  find  the  solution.  The  value  0.01  seems  to  be  a  good  tradeoff  be- 
tween time  and  accuracy.   Euler's  Method  can  be  stated  in  the  following  gen- 
eral form.  y    =  y  +  vy.^t 

n+1   ■'n   ^ 

To  solve  the  differential  equation  y  +  y  +  f (t)  where  f(t)  is  the  input,  the 

following  difference  equations  are  solved  repeatedly. 

t  ,  ,  =  t  +  At 
n+1    n 


y  =-y  +f(t) 
•'n   •'n    ^    ' 

n+1   •'n   •'n  ** 
The  method  can  be  generalized  to  solve  an  W^   order  differential  equation 
by  changing  it  to  a  set  of  N  first  order  differential  equations. 


CHAPTER  III 


PROCEDURE 


The  description  of  the  program  can  be  divided  into  six  main  parts. 
These  are  Input,  Phase  I,  Phase  2,  Phase  3,  Output,  and  System  Simulation. 
Through  the  input  the  user  can  specify  the  order  of  the  differential  equa- 
tion, the  initial  parameter  values  and  the  initial  stepsize  values.   By 
changing  the  appropriate  FUNCTION  Subroutines  the  user  can  specify  the  input 
f(t)  and  the  distance  function.   Phases  1  and  2  are  used  to  insure  the  con- 
vergence of  the  response  so  there  will  be  a  solution  for  the  settling  time. 
Phase  3  makes  the  final  parameter  adjustments  and  makes  adjustments  to  dis- 
continuities in  the  solution  space.   A  simplified  distance  function  defined 
as  the  integral  from  zero  to  twenty-five  seconds  of  the  absolute  value  be- 
tween the  input  and  output  is  used  in  Phases  1  and  2.   This  is  a  continuous 
function  that  simplifies  the  search  procedure.   Phase  1  uses  steepest  descent 
search  technique.   Phases  2  and  3  use  a  relaxation  search  method.   Each  time 
the  distance  function  is  evaluated  pertinent  information  is  printed  out  so 
the  user  can  know  where  the  search  procedure  has  been.   When  the  optimal  par- 
ameters are  found  they  are  printed  out  with  the  message  PARAMETERS  OPTIMIZED. 
The  system  simulation  obtains  the  necessary  information  for  calculating  the 
distance  function  and  uses  a  generalized  Euler's  Method  for  solving  N^''  order 
differential  equations. 

The  program  is  basically  written  in  FORTRAN  II.   The  single  exception  to 
this  is  the  IMPLICIT  statement  which  makes  all  floating  point  variables 


>-i^*".j>iH*'»'V"*'  ' 
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double  precision.   This  means  that  every  time  the  ABS,  or  SQRT  Library  Sub- 
routines are  used  they  must  be  specified  as  DABS  and  DSQT.   The  program  was 
run  on  an  IBM  360/50  using  the  FORTRAN  IV  Level  G  compiler  on  Release  12. 
The  execution  time  for  a  typical  run  on  the  third  order  problem  described  ' 
in  this  paper  is  6  minutes.  This  time  increases  with  increasing  order  of  the 
differential  equation  and  decreasing  stepsize. 

INPUT 

The  user  must  specify  the  order  of  the  differential  equation,  the  intial 
parameter  values,  and  the  initial  stepsizes.   The  program  has  a  default  so 
that  if  the  user  has  no  idea  what  values  the  parameters  should  be  the  program 
will  arbitrarily  set  them  to  1  and  the  initial  stepsize  to  one  tenth.   JJ  re- 
presents the  order  of  the  differential  equation  and  KJ  is  the  default  code. 
JJJ  is  set  equal  to  JJ-1.   This  is  the  number  of  adjustable  parameters.   lERR 
is  a  code  that  determines  which  ERROR  definition  is  to  be  used  (see  page  33). 
AA  is  the  time  of  integration  in  seconds  used  in  Phase  1  and  2.   DELTX  is  the 
time  increment  used  in  Phase  1  and  2.   The  input  list  JJ,KJ, IERR,AA,DELTX  is 
read  in  on  a  3ll,2X,2F5.0  format.   JJ  must  be  assigned  a  value  (2<JJ<8)  but 
all  the  other  variables  will  take  on  a  default  value  if  left  blank.   The  de- 
fault values  are  KJ=0,  IERR=1,  AA=25.,  DELTZ=.05.   If  JJ  is  9  the  program 
will  terminate.   If  KJ  is  positive  two  more  cards  are  read  with  a  8F10.5  for- 
mat.  The  first  card  should  have  parameter  values  A.  in  ascending  order.   The 
second  card  should  have  stepsizes  DELT.  in  ascending  order.   If  KJ  is  positive 
the  program  goes  directly  to  Phase  3.   The  default  occurs  if  KJ  is  zero  or 
blank.   The  input  f(t)  and  distance  furction  are  modified  by  changing  the  ap- 
propriate FUNCTION  Subroutines.   (see  Appendix  B) 
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PHASE  1 

This  Phase  used  steepest  descent  to  adjust  the  parameters  A.  from  their 
initial  value  of  1  to  a  set  of  parameters  whose  response  converges  to  the  in- 
put with  minimal  error.   All  parameters  are  adjusted  simultaneously  along  a 
straight  line  tangent  to  the  gradient  at  the  starting  point.   This  continues 
until  a  local  minimum  is  reached  or  the  parameters  are  adjusted  15  times. 
The  parameter  values  are  taken  as  a  new  starting  point  and  the  process  is  re- 
peated two  times.   DELT.  is  defined  to  be  the  stepsize.   The  interval  of  in- 
tegration, AA,  is  25  seconds.   When  the  DELT.  are  calculated  only  one  para- 
meter is  varied  at  a  time.   ERR  is  the  absolute  value  of  the  difference  be- 
tween the  input  and  output  ]y  -  f (t)]  . 


STEP  1 

The  following  quantities  are  evaluated. 


-  S: 


ERR  dt 


ERR  dt 


for  A.  -  DELT.  /lO 
1       1 

for  A. 


for  A.  +  DELT.  /lO 
1       1 


H3=   JeRR   dt 
From  these  values  four  cases  are  obtained  as  follows: 
The  resulting  curves  are  shown  in  Fig.  2. 
CASE 


HI 


H2 


lA 

Hl>  H2>  H3 

IB 

Hl<  H2<  H3 

IC 

H2>  HI,    H3 

ID 

H2<  HI,    H3 

H3 


Fig.  2  Sample  Curves  Showing  Cases 
lA,  IB,  IC,  and  ID. 
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PI  IAS K  1 
I 

N  =  0 
1  =  0 


T 


1  =  1  +  1 


CALCULATE 
AND  TEST 
HI   H2   H3 


CASE  1 


CASE  2 


CASE  3 


CASE  4 


DELT ( I ) =F/ABS (DELT ( I ) ) 
.005£DELT(I)  S  .38 


A(I)=  min(Hl,H3 


DELT(I)/10 


A  li  li 


no 


DIST(1)=H3 


P(I)=P(I)+DELT(I)   1=1, JJJ 


DIST(2)=  /err  dt 


yes 


DIST  =DIST2 


P(r)=P(l)-DELT(l)   1=1, JJJ 
N  =  N  +  1 


no 


RETURN 


Fig.  3   Flowchart  of  Phase  1 


F  is  defined  as  (ill  -  112). 

After  testing  to  sec  which  case  exists  the  following  action  is  taken. 

Case  lA  and  ID   DELT.  is  redefined  as  F/  DELT..   The  program  then  goes 
to  sfep  2.  ^ 

Case  IC         A.  is  redefined  to  that  parameter  value  corresponding 
to  the  smaller  of  HI  or  H2.   STEP  1  is  then  repeated. 

Case  ID         DELT.  is  decreased  by  an  order  of  magnitude  and  STEP  1 
is  repeated 


STEP  2 

The  resulting  DELT.  is  tested  as  follows: 

If    DELT.  >  0.38  DELT.  is  set  equal  to  0.38 

1       ■  1         ^ 

If    DELT.  <  0.005  DELT.  is  set  equal  to  0.005 

1  1         ^ 

If  all  the  DELT.  have  been  calculated  the  program  goes  to  STEP  3.   If 


not  I  is  incremented  by  1  and  the  program  goes  to  STEP  1. 


STEP  3 


After  all  the  DELT.  are  calculated  each  A.  is  incremented  by  the  corres- 
ponding  DELT..   ^ERR  dt  is  then  evaluated  repeatedly,  incrementing  the  par- 
ameters each  time.   When  the  value  of  the  integral  is  larger  than  the  pre- 
vious value  a  local  minimum  has  been  reached.   All  A.  are  then  decremented 
by  DELT^,  I  is  set  to  one,  and  the  counter  KJJ  is  incremented  by  one.   If  KJJ 
is  three  the  program  goes  to  Phase  2.   Otherwise  the  program  goes  to  STEP  1. 
See  Fig.  3  for  a  flow  diagram  of  Phase  1. 

PHASE  2 

This  phase  uses  a  relaxation  search  method  to  make  the  parameter  re- 
sponse converge  even  more  quickly.  A  single  parameter  is  varied  with  all 
others  held  fixed.   The  parameter  is  varied  until  one  of  three  things  occur. 


ii 


lA 


First,  the  dificrcncc  between  the  two  distance  functions  is  negligible,  se- 
cond, the  stepsize  may  get  too  small,  or  third,  the  parameter  has  been  varied 
more  than  15  times.   If  one  of  these  criterion  is  true  the  next  parameter  is 
then  varied.   Each  parameter  is  varied  in  turn  until  the  change  in  the  dis- 
tance function  for  JJJ  successive  parameter  changes  is  negligible.   The  dis- 
tance function  used  is  the  integral  from  zero  to  AA  seconds  of  the  absolute 
value  of  the  difference  between  the  input  f(t)  and  the  output.   AA  is  25  sec- 
onds in  this  report. 

The  initial  parameter  values  A.  are  those  calculated  in  Phase  One.   The 

initial  stepsize  DELT.  are  one  tenth.   All  changes  in  stepsize  are  by  orders 

^         1 

of  magnitude.  CODE  is  plus  or  minus  one  corresponding  to  positive  or  neg- 
ative slope  of  the  distance  function.   DIST(K)  corresponds  to  the  most  recent 
value  of  the  distance  function,  DIST(K-l)  the  previous  value,  etc.   Each  time 
the  simplified  distance  function  is  evaluated  it  is  done  twice,  once  for  A^, 

and  once  for  A  +  DELT."  CODE/10.   By  this  means  the  slope  can  be  determined 

i       1 

and  CODE  set  to  the  proper  value.   The  parameter  A.  is  then  incremented  by 
DELT."  CODE.   Two  more  distance  functions  are  then  evaluated  corresponding  to 

the  new  A  and  A  +  DELT."  CODE/IO.   When  this  is  done  one  of  the  following 

ill 

Cases  will  result.   CODE  is  taken  to  be  plus  one  in  this  example.   This  im- 
plies that  DIST(l)  is  greater  than  DIST(2).   See  Fig.  4  for  sample  curves. 

CASE   2A     DIST(1)>  DIST(2)>  DIST(3)>DIST(4) 

This  case  is  a  desirable  result.   It  means  the  search  is 
continuing  down  a  slope  in  the  parameter  space. 

CASE   2B     DIST(l) >DIST(3),   DIST(4) >  DIST(3) 

This  case  shows  the  minimum  has  been  crossed  and  that  it 
is  nearest  DIST(3). 


IP  jp'i.ii  'i'  ii-"»^-'_j^r'w  *^vi.'j«i  i.'.w«!..  I ' 
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CASE  2C     DIST(3)>  DIST(l),    DIST(4) >  DIST(3) 

This  case  shows  the  minimum  has  been  crossed  and  that  It  is 
probably  nearest  DIST(2). 

CASE   2D     DIST(3)>DIST(1),     DIST(3)  >  DIST(4) 

This  case  indicates  that  a  discontinuity  has  been  crossed 
in  the  solution  space  or  that  the  stepsize  is  very  large. 


DIST 


CASE  2A 


CASE  2B 


CASE  2C 


CASE  2D 


Fig.  4  Sample  Curves  Showing  Cases  2A,  2B ,  2C,  and  2D. 
The  above  cases  are  tested  for  and  the  following  action  is  taken. 
CASE   2A   A.  is  incremented  by  DELT  . 

CASE   2B   DELT^  is  decreased  by  an  order  of  magnitude  and  3*C0DE>'fDELT . 

is  subtracted  from  A. .  ^ 

1 

CASE   2C   DELT^  is  decreased  by  an  order  of  magnitude  and  7*C0DE>'fDELT. 

is  subtracted  from  A..  ^ 

1 


CASE   2D   DELT 
the  p 


^  is  decreased  by  an  order  of  magnitude  and  A.  is  set  to 
parameter  value  corresponding  to  DIST(2).    ^ 

If  one  of  the  following  criterion  are  not  satisfied,  the  program  sets 

DIST(1)=  DIST(3)   and  DIST(2)=  DIST(4).   The  program  then  sets  CODE  to  the 

proper  value  or  finds  DIST(3)  and  DIST(4)  corresponding  to  the  new  A  and  A 

i      i 

DELT(l)^vcODE/10. 

1.    |disT(3)-DIST(4)|  <DMIN    This  means  negligible  improvement  from 

the  change  in  parameter. 
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I        i'liASK    1 

T  =  L 

DALOW=.0001 

DTKKM=. 00005 

IW I N=. 000001 

C()D1';=1 

M=l 

ICONT=0 


IC0NT=IC0NT+1 


no 


CALL  C 

A .  =A .  +DELT .  -'.-CODE/10 
11  1 

CALL  C 

A .  =A .  -DELT .  -"-CODE/ 1 0 
11  1 

F=DIST    -DIST 

G=DIST^-DIST^ 


DKi.T  =r)F';L'iyvio 

LCONT=0 

PRINT 

DELT.=.l 

D.=DtsT„ 
1  3 

I=T-H 


I-JJ-1 

/^ 

1 

1  =  1 

-  \/(J 

M=I 

. 

JJJ 


M=2 


Fig.  5   Flowchart  of  Phase  2 


DELT   <DALOW 


3.    IC0NT>15 
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This  indicates  that  the  stcpsizc  is  less 
than  that  allowed  by  the  user. 

This  means  a  single  parameter  has  been 
adjusted  more  than  15  times  successively. 
Levine(1965)  shows  that  by  limiting 
the  number  of  adjustments  of  a  single 
parameter  the  minimum  is  found  more 
quickly. 

When  one  of  the  above  criterion  is  satisfied  D.  is  defined  to  be  equal 

1 

to  the  last  DIST.  calculated.   The  following  overall  termination  criterion 
is  then  tested. 


JJJ 

If    T~   Id.  .-d. 

h-y,     I  1-1   1 


i=2 


<  DTERM 


The  program  goes  to  Phase  3. 


If  this  inequality  is  not  satisfied  I  is  changed  so  the  next  parameter 
will  be  adjusted.   The  program  returns  to  set  CODE  to  the  proper  value  and 
the  process  is  repeated  until  the  termination  criteria  is  satisfied.   See 
Fig.  5  for  a  flow  diagram  of  Phase  2. 

PHASE  3 


This  phase  is  very  similar  to  Phase  2.   However  there  are  three  impor- 
tant differences.   First  the  required  difference  function  Jt  +ERROR   is 
used,  second  a  smaller  stepsize  DELTX  is  used,  and  third   some  additional 
tests  are  made  to  determine  when  the  search  has  crossed  a  discontinuity. 
When  a  discontinuity  is  detected,  the  stepsize  is  decreased  by  an  order  of 
magnitude  and  the  parameter  is  adjusted  so  it  is  on  the  small  side  of  the 
discontinuity.   If  the  termination  criterion  is  satisfied  while  on  the  wrong 
side  of  a  discontinuity,  A^  is  set  to  its  previous  value  and  D.  is  set  equal 
to  DIST^  instead  of  DIST  .   The  action  taken  after  testing  for  the  four 
cases  (see  Phase  2)  is  different  for  cases  3B  and  3C. 


'IIASK  3 
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1=1 

DALOW=.0001 

DTERM=. 00005 

DMIN=. 000001 

P1=0. 

P2=0. 

ML=0 

C0DE=1 

M=l 

ICONT=0 


CALL  C 

P1=P2 

P2=A. 

Ai=AJ;+DELT^'VCODE/l0 

CALL  C 

Ai=Ai-DELT^>vcODE/10 

F=DISTo-DIST, 


A.=A.+DELT.VcCODE 

1      1  1 


A.  =A.  -DELTi-'.-CODE'VlO 
1       ^ J: 


M=l 


Ai=Aj^+DELTj_-A-CODE 


A£=A^-DELTi*C0DEVf2 

DELTi=DELTi/10 

M=l 


Fig.  6  Flowchart  of  Phase  3 

continued  on  next  page 


19 


DF,LT.=DELT.*10 
PR  I  N't     ^ 
ICONT=0 


I-JJ+1     >— 

\y^ 

1=1 

—     1^  0 

1 

1 

M=1 

Fig.  6  Continued  Flowchart  of  Phase  3 
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Case   3B   DELT.  is  decreased  by  an  order  of  magnitud 


e. 


Case   3C   DELT,  is  decreased  by  an  order  of  magnitude  and 

10Vcc4dEV.-DELT.  is  subtracted  from  A  . 
1  1 

When  Phase  3  is  terminated  the  program  prints  PARAMETERS  OPTIMIZED 
and  prints  out  the  optimal  parameter  values.   See  Fig.  6  for  a  flow  diagram 
of  Phase  3. 

OUTPUT 

The  program  will  terminate  in  one  of  three  ways.   The  desired  termina- 
tion is  taken  when  the  parameters  are  optimized.   If  the  parameters  have  a 

response  that  does  not  converge  within  25  seconds  the  program  will  also  term- 

I  2      2" 
mate.   It  the  distance  function  -^T  +ERROR   is  evaluated  more  than  400 

times  for  the  same  input  data,  the  program  will  terminate  and  print  out  the 
message  PROGRAM  TERMINATED  AFTER  401  LOOPS.   Because  of  the  possibility  of 
these  undesired  terminations   the  user  needs  to  know  what  the  program  has 
done  up  to  the  time  of  the  termination.   Each  phase  prints  pertinent  infor- 
mation applicable  to  that  Phase. 

In  Phase  1  the  new  DELT^  are  printed  each  time  they  are  evaluated. 
Each  time  the  parameters  are  incremented  the  program  prints  out  the  corres- 
ponding distance  function  and  the  parameter  values.   In  Phase  2  the  program 
prints  DIST3,  DIST^,  CODE,  I,  A.,  and  DELT..   In  Phase  3  the  program  prints 
the  same  as  in  Phase  2  plus  the  values  for  AREA  and  T.   Whenever  I  is  incre- 
mented in  Phases  2  and  3,  the  parameter  values  and  DELT.  are  printed.   At 
the  end  of  the  program  if  the  proper  termination  criterion  is  satisfied  the 
program  prints  PARAMETERS  OPTIMIZED  and  prints  the  parameter  values. 

SYSTEM  SIMULATION 
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The  actual  mathematics  of  the  simulation  is  described  in  the  theoretic- 
al analysis.   Interwoven  with  this  are  all  the  necessary  checks  to  evaluate 
the  distance  functions.   Some  problems  are  encountered  here  because  of  round 
off  error  and  the  inability  of  a  binary  computer  to  represent  decimal  numbers 
exactly.   Double  precision  is  used  on  all  floating  point  variables  to  try  and 
improve  the  accuracy.   The  round  off  error  can  be  most  readily  seen  in  the 
value  representing  time.   Because  it  is  determined  by  repeatedly  adding  small 
increments  a  large  number  of  times  the  roundoff  error  becomes  very  apparent. 
This  also  affects  the  accuracy  of  the  solution  to  the  differential  equation. 

If  the  value  of  time  exceeds  30  seconds  when  using  the  required  distance 
function,  T   is  set  equal  to  25,  the  distance  function  calculated,  and  control 
is  returned  to  Phase  3. 

LLL  is  the  key  used  to  determine  which  distance  function  is  to  be  eval- 
uated.  LLL  is  equal  to  one  for  Phases  1  and  2.   This  causes  the  subroutine 
to  evaluate  the  continuous  distance  function.   LLL  is  equal  to  two  in  Phase  3. 
This  causes  the  assumed  distance  function  to  be  used.   LLL  is  equal  to  three 
when  a  listing  is  made  of  the  response  of  the  differential  equation  versus 
time. 


CHAPTER  IV 


RESULTS  AND  CONCLUSION 


RESULTS 

The  problem  to  be  solved  is  repeated  as  follows.  The  parameters  of  a 
third  order  differential  equation  are  to  be  found  which  minimize  the  distance 
function  -aT  +ERROR  with  a  modified  step  input.   T   is  the  time  required  for 
the  response  to  decrease  to  and  stay  within  five  percent  of  the  final  value 
of  the  input.   ERR  is  defined  to  be  the  difference  between  the  input  and  the 
output.   Three  different  ERROR  definitions  are  used  as  follows. 

1.  lAE   Integral  of  the  absolute  value  of  ERR. 

2.  lES   Integral  of  ERR  squared. 

3.  ITES   Integral  of  time  times  ERR  squared. 

In  addition  to  solving  the  above  problem,  optimal  parameter  values  were 
found  for  the  simplified  distance  function  ERROR  for  each  of  the  three  defi- 
nitions.  This  was  done  for  a  comparison  of  the  performance  measures  of  the 
distance  functions  for  each  ERROR  definition.   See  Table  1  for  the  tabulation 
of  the  performance  data  and  optimal  parameter  values. 


For  clarity  the  distance  functions  -JT^+ERROR^  and  ERROR  are  hereafter 
referred  to  as  DFl  and  DF2  respectively. 

Because  of  the  definition  of  DFl  the  optimal  parameters  are  the  same  for 
all  three  ERROR  definitions.   This  result  is  caused  by  the  high  sensitivity 
of  DFl  to  the  settling  time  T^  for  this  order  differential  equation.   There 
is  only  one  set  of  parameters  that  cause  the  output  to  respond  most  quickly 
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with  no  more  than  five  percent  overshoot  or  undershoot.   Because  of  this  the 
values  of  ERROR  have  little  or  no  effect  in  determinimg  the  optimum  parameter 
values  with  DFl.   See  Fig.  7  for  a  plot  of  the  optimal  response  versus  time 
for  DFl. 

When  DF2  is  used  there  are  no  constraints  on  overshoot  and  undershoot. 
The  time  of  integration  must  be  specified  as  twenty-five  seconds  because  it 
is  not  implied  in  the  definition  of  DF2  as  it  is  in  DFl.   The  performance 
values  for  DF2-IAE  are  nearly  the  same  as  those  obtained  from  DFl.   DF2-IES 
and  DF2-ITES  have  faster  response  times  but  the  settling  time  and  overshoot 
are  increased.   From  the  data  in  Table  1  it  can  be  seen  that  when  one  per- 
formance measure  improves  one  or  more  of  the  others  get  worse.   Thus  there 
is  a  tradeoff  between  performance  measures  and  good  engineering  judgment  must 
be  used  to  select  the  distance  function  that  emphasizes  the  most  desirable 
combination  of  performance  measures.   See  Fig.  8  for  a  plot  of  the  optimal 
response  versus  time  for  the  three  ERROR  definitions  for  DF2. 

A  small  change  in  the  stepsize  DELTX  can  cause  rather  significant 
changes  in  the  calculated  optimum  parameter  values.   This  is  caused  by  the 
accuracy  limitations  of  Euler's  Method  with  large  stepsize  values.   If  higher 
accuracy  is  needed  some  other  method  such  as  Runge-Kutta  should  be  used  for 
solving  the  differential  equation.   In  most  cases  however,  high  accuracy  is 
not  needed.   With  good  judgment  fairly  large  stepsizes  can  be  used  to  deter- 
min  rough  optimum  values  with  very  little  computation  time. 

The  weakest  part  of  the  program  is  the  time  required  to  solve  the  diff- 
erential equation.   Since  the  system  must  be  simulated  many  times  in  order  to 
find  the  optimum  values,  most  of  the  computer  time  is  used  in  doing  this  sim- 
ulation.  In  order  to  reduce  the  cost  of  finding  the  solution  the  simulation 


Ik 


01 


on 

< 
w 
S 

w 
u 

^ 

pi 
o 

w 
a. 

< 

u 
> 
w 
to 

cc 
o 

CO 


(is 

w 

K 

kJ 

P 

m 

<; 

Qi 

H 

o 

pi 

lA 

M 

to 
o 

M 
O 


Ixl 
O 


O 

z 
o 
CO 


o 


pi  w 
O   3 


Pi 
Pi 

w 


o 


S 


2: 


CO 

J  g   o 

HMO 
W  to 

to 


CO 

to  S    O 

M    h- 1      O 

Pi  H    <u 

CO 


w  <;  d 

2  hJ  o 
M  w  o 
9-1  p    01 

CO 


O 
O 

to 

a: 
u 
Q 
Z 
3 


O 
O 

X 

to 

Pi 

> 
o 


4-) 

c 
a) 
o 
u 

d) 


4-1 

C 

(U 

o 

(-1 

a, 


CO 

Pi 
w 

w 

s 


o 

w 

Pi 

(X 

Pi 

>- 

w 

H 

w 

Z 

C) 

O 

z 

h-l 

< 

H 

H 

C) 

CO 

Z 

M 

3 

G 

fc 

C7>  r-^  o^  vO  O 

<r        o        c^       "~i        o 

CN  O'         O  --<  c^ 


vD 


CX) 


o       o       o 

CM  (M  CM  CS 


.-<  CX) 


CM 

0^ 


en 
CM 


CO 
CN 


CM 


in 


CO 


CX) 

ON 


00 


<X) 


CO 

in 


o 

CX) 


o 

o 


o 

o 


o 
o 


o 


1^ 


in 


in 


in 
t— ^ 
O 

in 
o 

in 

r-l 

o 

o 

CO 

o 

00 

ON 

CM 

ON 
ON 

CM 

Cvl 

CM 

CM 

I— f 

I— 1 

• 

• 

o 

• 

o 

• 

ON 
ON 

• 

00 

CO 

CO 

w 

to 

w 

w 

CO 

w 

< 

w 

H 

< 

w 

H 

M 

M 

M 

M 

M 

M 

.— 1 

r-t 

CM 

CM 

CM 

Cn 

fa 

fe 

fc. 

fa 

p 

P 

P 

P 

P 

g 

O 

M 

•4-1 

• 

•o 

(U 

4-> 

CO 

U      ' 

W)   CO 

CU     TD 

4-j    d 

c  o 

•H      O 

<u 

CO    CO 

•H 

C 

CD     (U 

O.   (U 

X  ■!-> 

4-1     IJ-I 

•H 

M    <4-4 

o 

!-t     O 

M     4J 

CU 

O 

4)     >-i 

x:   CU 

CN 

H     N 

Pi 

O 

Pi 

pi 

w 

pi 

+ 

o 

CM     CO 

Pi 

H 

Pi 

W 

II 

II 

i-( 

CM 

fa 

fa 

P 

P 

25 


picxjcocuozcow 


26 


ii 


o 

o^ 

<NI 

CO 

00 

ON 

o 

ON 

o 

\ 


o 


00 

t— I 

O  .-I 


0^ 
ON 


3 

w 

a. 

w 

CO 

w 

C 

< 

w 

fc- 

0^ 


00 


CO 

C 

o 
o 

to 


CM 

o 

M 
O 


I 

•H 

CO 

D 
CO 
u 

> 


a, 
P 


« 

(U 
(0 

c 
o 

Q. 
W 

(U 

OS 


o 
•-t 

CM 

00 

00 
•i-l 


eiiwcocLj02c«w 


27 


routine  must  either  be  used  less  or  a  cheaper  method  of  simulation  used. 
Better  prediction  or  extrapolation  methods  would  help  reduce  the  number  of 
times  the  routine  must  be  used.  A  high  order  system  could  be  simulated  more 
cheaply  on  a  hybrid  computer.   For  this  reason,  the  program  was  written  so 
it  could  be  adapted  to  a  hybrid  computer. 

CONCLUSION 

The  program  presented  in  this  paper  is  a  practical  approach  to  the 
optimal  parameter  design  problem.   With  the  use  of  the  program  and  a  judi- 
cious choice  of  a  distance  function  that  emphasizes  the  desired  performance 
measures,  the  user  can  quickly  determine  the  optimal  parameter  values  for 
an  N^*^  order  (2i  N<8)  differential  equation.   The  distance  function  should, 
if  at  all  possible,  be  chosen  so  it  is  a  continuous  function.   Only  when 
the  system  must  have  special  response  characteristics  should  a  discontinuous 
distance  function  be  used. 

The  next  step  in  the  development  of  the  program  should  be  to  generalize 
it  so  the  program  could  simulate  systems  with  both  poles  and  zeros.   This 
would  greatly  increase  the  usability  of  the  program. 
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APPENDIX  A 
DEFINITION  OF  PERFORMANCE  MEASURES 


These  terms  are  defined  somewhat  differently  than  usual  because  they 
are  not  limited  to  step  inputs.   Any  input  f(t)  can  be  used  that  has  a  con- 
stant final  value  after  some  time  T-,  and  never  exceeds  the  final  value. 


Overshoot  Percent  overshoot  is  defined  as  the  difference  between 
the  maximum  response  r(t)  after  T-^  and  the  final  value 
of  f(t)  times  100  divided  by  the  final  value. 


Undershoot       Percent  undershoot  is  defined  as  the  difference  between 
the  final  value  of  f(t)  and  the  minimum  response  after 
the  first  negative  slope  crossing  of  the  final  value 
line  after  T, ,  multiplied  times  100  and  divided  by  the 
final  value. 


Time  Delay       Time  delay  is  defined  as  that  time  Tj  required  for  the 
response  to  reach  50  percent  of  its  final  value. 


Rise  Time        The  rise  time  T^  is  defined  as  the  time  required  for 
the  response  to  rise  from  10  percent  of  its  final 
value  to  90  percent  of  its  final  value. 


Settling  Time    The  settling  time  Tg  is  defined  as  the  time  required 
for  the  response  to  decrease  to  and  stay  within  5 
percent  of  its  final  value. 
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APPENDIX  B 
COMPUTER  PROGRAM 

The  following  pages  contain  a  listing  of  the  Fortran  program  described 
in  this  paper.  The  program  consists  of  a  main  calling  routine,  four  subrou- 
tines, and  three  FUNCTION  subroutines.   The  program  was  written  in  this 
manner  so  the  user  could  easily  change  it  to  meet  his  needs.   If  a  faster 
execution  time  is  desired,  the  FUNCTION  subroutines  should  be  incorporated 
into  the  calling  routine.   Comment  statements  are  interspersed  throughout 
the  program  to  help  the  user  understand  what  the  program  does. 

The  following  are  explanations  of  the  three  FUNCTION  subroutines. 

FUNCTION  AINPUT(Z) 

This  routine  calculates  the  value  of  the  input  as  a  function  of 
time.   Z  is  the  variable  that  represents  time.   The  user  can  change 
this  routine  as  desired. 

FUNCTION  DISTFN(AREA,T) 

This  routine  calculates  the  value  of  the  distance  function  called 
for  in  the  problem  described  in  this  paper.   If  this  routine  is  changed 
to  use  other  values  in  the  distance  function,  SUBROUTINE  C  should  be 
changed  accordingly  to  supply  needed  values  or  delete  uneeded  values. 
The  user  can  change  this  routine  as  desired. 

FUNCTION  RRR(ERR,Z) 

This  routine  is  used  to  allow  several  different  ERROR  definitions 
to  be  evaluated  in  the  same  program  run.   The  program  keys  off  the  value 
of  lERR  to  determine  the  ERROR  definition.   Care  shoula  be  taken  when 
changing  this  routine  because  of  the  COMMON  statement. 


//CPTIMI2F  JC  i       0?F/t0A08G00?»-3WlTZER-»MSGLEVFL  =  l       OlO»OAO»000 

//  EXEC  FTGCLGKS.PARM.FCRT=-BCD- 

//FCRT.SYSIN    DD    ♦ 

C 

C 

C  ***********************  4;  *-^f  *  ^( -1^ ^--K  ^{■»^^^;•  ■«**«***** *******«*»«»#i( *#»*♦**#♦ 

C  ***PARAMET€R  OPTIMIZATION  OF  A  JJ  TH  ORDER  DIFFERENTIAL  EQUATION** 

C  ***KFNNETH  SWITZFR    1968*** 

C  **«****-^**«K  ***«•*  ^<•*****^<'*•***^;*-^^«*»***•»*^f***^t*^<.^^s.)(.^^.^^»^^^(^^^^#^(^^^^^f^^^^>^^^^f 

r  **      THIS  PROr-RAM  WILL  FIND  THF  OPTIMUM  PaRAMFTFRS  SURJFCT  TO  ** 

C  **  A  SPECIFIED  DISTANCE  FUNCTION  AND  INPUT.   THE  INPUT  MUST  HAVE** 

C  **-  A  CONSTANT  FINAL  VALUE  AFTER  SOME  TIME  TO  AND  NEVER  EXCEED  ** 

C  *^^  THAT  FINAL  VALUF.   THF  DISTANCE  FUNCTION  IS  THE  SQUARE  ROOT   ** 

C  **  OF  THE  SUM  OF  THE  SQUARES  OF  THE  SETTLING  TIME  T(S)   AND  THE  ** 

C  **  VALUE  OF  THE  INTEGRAL  FROM  ZERO  TO  T(S)  OF  SOME  FORM  OF  THE   ** 

C  «■*  ERROR  BETWEEN  THE  INPUT  AND  THE  OUTPUT.   T(S)  IS  DEFINED  AS   ** 

C  **  THE  TIME  REQUIRED  FOR  THF  RFSPONSE  TO  DECREASE  TO  AND  STAY    ** 

C  **  WITHIN  .C5   OF  THE  FINAL  VALUE.   IF  lERR  IS  1  EER  IS  EQUAL  TO** 

C  **  ERR.  IF  lERR  IS  ?  FFR  IS  ERR  SQUARED.   IF  IFRR  IS  3  EER  IS    ** 

C  **  T  TIMES  (^RR  SQUARED.   THF  APPROPRIATE  FUNCTION  SUBROUTINES    ** 

C  **  CAN  RF  CHANGE  TO  SUIT  THF  USERS  NEEDS.   EULFRS  METHOD  IS  USFD** 

C  ^<    TO  SOLVE  THE  DIFFERENTIAL  EQUATION.   THE  PROGRAM  GOES  THROUGH** 

C  **  TWO  PHASES  TO  INSURE  A  SOLUTION  TO  T(S)  IN  PHAS^  3.  ** 

C  **       INPUT  #* 

r  **  THE  ONLY  VALUE  REQUIRED  IS  JJ.   ALL  OTHER  VALUES  WILL  DEFAULT** 

C  **  THE  INPUT  IS  OF  THE  FORM  J J  ,  K J  ,  I  ERR , AA  ♦  DELTX .   THIS  IS  READ   ** 

C  **  ON  A  3n,2F5.0  FORMAT.   IF  JJ  IS  9  THE  PROGRAM  WILL  TERMINATE** 

C  ^^  DEFAULT  VALUES  IF  THE  LOCATIONS  ARE  LEFT  BLANK  ARE  AS  FOLLOWS** 

C  **  KJ  =  l>    IERR  =  1     AA  =  ??.    DFLTX  =  .OS.    KJ=1  IS  A  CODE  USED  IF** 

C  **  THF  STARTING  POINT  IS  TO  RE  SPECIFIED.   TWO  MORF  CARDS  ARE    ** 

C  **  READ  IN  WITH  A  aFlO.O  FORMAT.   THE  FIRST  CARD  MUST  CONTAIN    ** 

c  **  iH^   Parameter  values  a(i)  in  ascending  order,  the  second    ** 

C  *^  CARD  MUST  CONTAIN  THE  INITIAL  STEPSIZE  IN  A  SIMILAR  MANNER.   ** 

r  ******■^****^<^*-"-»-'«^^■*>i***^^*^^^^^<^^^t^**K«»^^*♦*-s.^(*Jt^(.^^^^(^(.^f^^(.^f^f^^^{.^J.^^^^^^^^^^ 

C 

^  *^'  PROGRAM  VARIABLES  *# 

^  "*  AA       INTEGRATION  TIME       DELT(I)  STEPSIZE  ** 

C  **  JJ       ORDER  OF  DIFF  EQ       KJ       CODE  FOR  DATA  READ  IN   ** 

C  ■<■*    Ml)  PARAMETERS  OF  DIFFERENTIAL  EQUATION  ** 

C  **  ALOW      vLOWARLF  DEVIATION  FROM  INPUT  FOR  SETTLING  TIME       ** 

C  **  IFRR     CODE  FOR  RRR  FUNCTION  SUBROUTINE  ** 

C  **  LLL      CODE  FOR  WHICH  PHASE  IS  USING  SUBROUTINE  C  ** 

C  **  LKK      CODE  FOR  THF  STATUS  OF  SETTLING  TIME  ** 

C  **  KKK      rCUNTPR  FOR  MAXIMUM  NUMBER  OF  DISTANCE  FN  EVALUATIONS** 

C  **  DFLTX    TIME  INCREMENT   THIS  IS  REDEFINED  LATER  FOR  PHASE  3   ** 

C  **  DIST(I)  THE  FOUR  STORED  DISTANCE  FUNCTION  POINTS  ** 

C  **  D(I)     USED  IN  TERMINATION  OF  PHASES  2   AND  3  ** 

r  ******^*****  ******  ******i;«****^)tic<r-^*>.(}^.**^-;t^Hj.^i.^t^^f^^t^(^f^^f^^^^^^^^^^^ 
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«*       MAIN  CR  CALLING  PROGRAM 
IMPLICIT  REAL^'aCA-H.vO-Z) 

niMrNf. iCN  Y(  i>) )  ,n(  1  c) )  ,A  ( ]() )  »nELT  ( 1  u)  ,niM  ( 0 ) 

COMMON  AA»ALOW,nrLTX,n»nF»  IFRR,LKK,,LLL,J,JJ»JJJ,<<<,K,  I 

**    INPUT  ARFA   AND  DFFAULT  CHFCK 
SO    Rf An98,JJ,KJ,IFRR»AA»DFLTX 

PRINTIC5 

PRINT  99,JJ,IFRR  ■    ■   " 

QO  F0RMAT(?I5') 
9i     FORMATOn  .2X,?F5.0) 

IF( JJ)80,£  D,81 
81  IF(JJ-9)83. 82.82 
S?    STOP 

83  IF(  IFRR)84.84»85 
64  IFRR=1 
8  5  IF(AA)88,88,89 

88  AA=25. 

89  IF(DFLTX) 92,92,87 
^2  DFLTX=.05 

87  no  809  1  =  1  ,JJ 

A  (  I  )=3. 

DFLT( I )=.l 
809  D( I )=0. 

K  =  U 

AL0W=.05 

LLL  =  1 

LKK  =  1 

<KK  =  0 
105  FORMAT( IHl ) 

no  743  I  1=]  ,4 
74^  nT?T(  I  n  =  ". 

A { JJ)=1. 

IF(KJ)802,802,803 
80-^  RFADIOO,  (A(  I  )  ,  1  =  1  ,JJ) 

RFADlOv.,(DFLT(  I  )  ,1  =  ]  ,JJ) 
ICO  F0RMAT(8F1C,5) 

GO  TO  1 
802  J=JJ+1 

JJJ=JJ-1 

CALL  PHASl ( A,DIST,DELT) 

DO  8011  =  1, J 
801  DFLT( I )=.l 

CALL  PHAS2 ( A,DIST,nrLT) 

PRINT  300  ■  . 

300  FORMAT(/,2X,31HSWITCH  TO  NFW  DISTANCE  FUNCTION,//) 

**  RFSFT  DFLTX  FOR  PHASE  3 

DFLTX=.G1 
]  LLL=2 

CALL  PHAS3{ A,DIST,DELT) 

LLL  =  3 

DE  =  0. 

**  CALL  FOR  LISTING  OF  RESPONSE  VALUES 
CALL  C( A,DTST,OFLT) 
GO  TO  8  0 
END 


#» 


#» 


♦« 


** 
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FLiNCTION    AINPUT(Z) 

IMPLiriT  RFAL*B( A-H,C-Z ) 
C      X*   THIS  FUNCTION  SUPRCUTINF  SPFCIFIFS  T^iF  INPUT  AG  A  FUNCTION   «« 
C      **     OF  Z  WHICH  HFPRESFNTS  TIME  IN  SECONDS  ** 

IF(Z-1.  )  1»2»2 
]  AINPUT=Z 

RFTURN 
?  AINPUT=1. 

RFTURN 

END 


FUNCTION  DISTFN(ARFA»T) 

J.VPLICIT  RFAL*8(A-H,0-Z) 
C      **  THIS  FUNCTION  SUBROUTINE  SPECIFIES  THE  DISTANCE  FUNCTION  USED** 
C      *•■*     IN  PHASE  3  IN  TERMS  OF  SETTLING  TIME  T  AND  AREA  *» 

DlSTf^N  =  DSQRT(AREA*ARFA  +  T*T) 

RFTURN 

END 

FUNCTION  RRR(ERR,Z) 

IMPLICIT  REALV.-8(A-H,0-Z) 

DIMENSION  Y(10),D(10),A(10),DELT(10),DIST(5) 

COMMON  AA,ALOW»OFLTX,r>,nF»IFRR,LKi<,LLL»J,JJ»JJJ»KK<,K,I 
C      **  THIS  FUNCTION  SUBROUTINE  KEYS  OFF  lERR  AND  LLL  TO  SPECIFY  THE** 
C      **      VALUE  OF  RRR  AND  IN  TURN  EER.  ♦♦ 

C      *•*  INTEGRAL  OF  ABSOLUTE  VALUE  OF  ERROR     lAE  «* 

C      ^-*    INTEGRAL  OF  ERROR  SQUARED  lES 

C      **  INTEGRAL  OF  TIME  MULTIPLIED  ERROR  SQUARED     ITES 

GO    TO     (  1  .2, 2), LLL 
1    RRR=FRR 

r^rjURN 
?    30    TO     C.A.'i)  »  IFRR 

•^  PPR  =  FRR 

RFTUR^' 

4  RRR=rRR*rqR 

RFTURN 

5  RRR=ERR*ERR*Z 
RETURN 

END 
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SUBROUTINF  PHAS ]  ( A  ,  D  I  ST  ,  OFLT )  • 

IMPLICIT  RFAL*8 ( A-H»C-Z ) 

nI^'FN?IC^J  A(  If.  )  ,COnF(]0)  ,0(10)  ,nFLl  (]0),rM5T{^),Y(lO) 

rOMKCN  AA,ALCW,orLTX,n»nr,TFRR,LK:<,LLL,J,JJ,JJJ,<.KK,K,I 
PHASf  1  USFS  A  STFFPf^^T  DFSCFNT  SEARCH  JFCHNIO'JF 
IINT=^RS     KJJ   MAXIMUM  ITFRATICNS  ALONG  LINF   KJ 


H2  H3 


C      ** 

c    **  c 

900  KJJ=0 
<J  =  0 
1=1 
C      **  EVALUATION  OF  HI 
1  CALL  C( A,niST,nFLT) 
CC=DABS(DFLT( I ) )/10. 
A ( I )=A{ I) +CC 
CALL  C( A,DIST,nFLT) 
A ( I )=A ( I) -p,vcr 
CALL  r( A,niST,PFLT) 

A ( I )-A ( n  +rc 

H1=DIST{K> 

H3=DIST(K-1 ) 

H?=DIST (K-2) 

F=H1-H2 

G=H2-H3 

H  =  F-G 
99  FORMAT( 15 ) 

PPIMT99, I 

PPINT71 tA { n 

PRINT??, HI ,H?,H^ 
71  FC^VATI  5H  A(  I)  ,F9.6  ) 
7?  FrRNiAT(  1X,?HH1  ,F12.8 

I F( F)40,OP ,5  0 
'^'O  IF(C)51  ,98,55 
51  '>FLT(  I  )=DFLT(  I)/10. 

GO  TO  1 
40  IF(G)55 ,98 ,42 
42  IF(Hl-i-t3)  47,47,48 

47  A{ I )=A( I ) +CC 
GO  TO  1 

48  A ( I )=A { I ) -CC 
GO  TO  1 

C      *^;  f^VALUATION  OF  nFLT(I) 
5*^  OFLTC  n  =F/nAPS(r^c-LT(  I  )  ) 

PPINT80,oc|_T(  I  ) 
C      **  CONSTRAIN  nFLT{ I ) 

IF(DABS(DE  LT(  I  )  )-.38)61  ,61  ,60 

60  DFLT( I )=.38*(F/DARS( F) ) 
GO  TO  300 

61  IF(DABS(DFLT( I))-. 005)62, 300, 300 

62  DFLT( I)  =  .0C5*(F/DABS(F)  ) 
300  PRINT80,DELT(  I  ) 

1  =  1  +  1 


MAX  RUNS 


«•* 


,?X,?HH?,F1?.8,2X,2HH3,F]2.8) 


?{ 

172 
171 

B 
80 


75 


81 

Oft 


IF( I-JJ) i ,70.70 

1  =  1 

**  INCREML  ^T  ON  STRAIGHT  LINE 

no  17]  I  T  =  l  ,JJJ 

A  {  I  n=A(  IT  )+OPLT(  I  I) 

CALL  C{ A,OIST,nELT) 

PRTNTRO,niST(K),(A(II),II=l,JJ) 

FORMAT ( 8F] 5.8 ) 

KJJ=KJJ+1 

IF(KJJ-15 )?8A,?73,?73 

I F{ DIST( K )-DrST(<-l  ))  172  ♦273,273 

KJJ  =  0 

DC  75  I  1  =  1  ,JJJ 

A  (  I  I  )=A(  II  )-OELT(  I  I ) 

**    TFRf^INATICN    CCHNTER 

KJ=KJ+1 

IF(KJ-3) 1,7A,74 

PPI  MT81 

-CRMAT{2X,6HSWITCH) 

RFTURM 

STOP 

END 
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SUF-ROUTINF  PHAf-?  (  A  ,  n  I  ST  ♦  n>FLT  ) 

IwPLiriT  Rf  AL.;<-8(  A-H,C-7  ) 

DIN' FN  SI  ON  niSTC^)  »nFL  T(  10)  ,A  (  10)  »n(  ]n) 

CCMMCN  AA,ALCW«DFLTX.D,DF,  IFRR,LKi<s,LLL,J.JJ»JJJ»fCKK.  ,K,I 
C      **    PHASF  7    USES  RELAXATION  FOR  OPTIMIZING  THE  PARAMETERS 
C      *<^   SEE  PHASE  THREE  FOR  AN  EXPLANATION  OF  SOME  OF  THE  VARIABLES 

1  =  1 

DALOW=.00  0l 

DTFPM=. 00005 

nviN!=.cooooi 

CCDE=i . 

y  =  1 

irONT=0 
1  irONT=irOMT+l 

IF(  irOMT-16)  3,-^0,-^0 
3  CALL  r(  A,r)IST,nFLT  ) 

A ( I )=A( I) +DElT( I)*C0DE/10. 

CALL  C(A,DIST»nELT) 

A( I)=A( I )-DELT( I )*C0DE/l0. 

-0RMAT(2X,7(F16.Q>?X)  j'SHCODF  ♦F3.n'»I5»2X»F16.9»2X,F16.0) 

PR  I  NT  1-7, 01  ST (K-1  )  »DIST(K)  ♦CODE 'I » A (  I  )  ,  DFLT ( I  ) 

GO  TO  (  10  ,70  )  ,M 

TF(niST(K-l  )-niST(K)  )1  1  ,1?,1 ? 

crnF=conp*(-i . ) 

IF(nAqS(niST(K-i  )-niST(K))-0MlM)30.1'.  ,1? 

A  ( i)=A  (  I)  +nFLT(  I  )*ror>F 

GO  TO  1        ■ 

IF(DIST (K-l)-OIST(K)  )21  »26»26 
DFLT( I )=DELT{ I) /lO. 
IF(DIST(;<-3)-DIST(K-l  )  )  2  3.  2  2,  22 

22  A  (  I  )=A( I  ) -  3ELT(  I)*C0DE*3.  •  • 
GO  TO  2A 

23  A (  I)=A( I)-DELT( I  )*rODF*7. 
2  4  M  =  l 

IF(DflT  (  I  )-r>ALOW)3n,3  0,7  5 
2  "5  IF{DIST(K  )-DIcT(K-1  )-DVIN)  3  0,30,] 
26  A  (  I  )=A(  I) +DELT(  n^CODF 

IF(nFLT{  I  )-D  ALOW)  •^0,2  7,27 
2  7  IF(DIST(.K-1  )-DIST(K)-DMIN)  30,30,1 


** 
** 


102 

06 

10 

1  1 

1  2 

1  ^ 

2C 

21 

37 


-^n 


10] 


40 
4] 

2 


A? 

43 

91 

9: 


PFLT( I )=nrLT{ I )*10. 

irc\'T  =  o 

PRINT!  ui  ♦n^^LK  I  )  ♦  (,A(  I  n  ♦11=1  ,JJ) 

F0PMAT(?X,F16.6»6F12.6) 

OFLT(  I  )=.] 

D( I  )=niST(K-l  ) 

1  =  1  +  1 

IF( I-JJ+1 )41 ,41 ,40 

1=1 

v  =  l 

DO    2    11  =  1, JJ 

DFLT(  I  I  )  =  ,]. 

SUN'  =  0. 

nn    4->    TI=?,JJJ 

SUM  =  ,SUN^  +  DAPS(n(  II-l  )-D(  II  )  ) 

IF(Suy-DT^Ry)43,l,l 

RFTURiN! 

ML=ML+i 

IF(ML-9) 1 ,92,92 

DFLT( I )=DFLT( I )*10. 

ML  =  C 

GO    TC    1 

FND 
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MAXIMUM  NUMBER  CF  INCRFMENTS 


;F  A(  I 


SUBROUTINE  PHAS3 ( A »D I  ST ,nELT ) 

IMPLICIT  RFAL*8(A-H»G-Z) 

niMFMSTCN  niST(  5)  »nFLT  (  10)  ,A  (  ]  0)  »[)(  10  ) 

rCMMCN  AA,ALCW»nFLTX,n,DF»IFRR»lXfC,LLL»J,JJ»JJJ»KKK,<,I 
C      «•*    PHAS!^  3  USFS  A  MCDIFFP  RELAXATION  METHOD 

1=1 
C      **  DALOW    ALLOWABLE  MINIMUM  STEPSIZE 
C      **  DTERM  PHASE  TFRMINATION  VALUE 
C      **  HMIN  PARAMETER  TERMINATION  VALUE 

nALOW=.OOC  1 

DTERM=.CGCG5 

DVIN=.UC0  0C1 

P  1  =  o  . 

P?  =  0. 

ML  =  C 

conF=i. 

M  =  l 
C      *-*  COUNTER  FOR 
ICONT=0 

1  IC0NT=IC0NT+1 

IE( IC0NT-16)3»30»30 
3  CALL  C( A,DIST,DELT) 
C      **  STORE  PREVIOUS  VALUE 

P1  =  P2 

P2=A( I ) 

A  (  I  )=A( I ) +nELT( I  )*C0DE/10, 

CALL  C( A,DlST,nELT) 

A(I)=A(I)-nELT(I )*C0DE/10. 

F0RMAT(7X,?(F16.9»?X ) ,SHCOnF  ,F3.G»I^,2X,E16.9,2X,F16.9) 

PRTNTl07,nrST(K-l  )  ♦r)IST(i<)  ,CODE»I  »A(  I  )  ,DELT(  I  ) 

GO  TO  ( 10,20)  ,M 

**  TEST  FOP  SLOPE 

IF(niST(K-l)-niST(K)  ) 11  ♦12^1? 

CODr=CODE*(-l. ) 

IE(DARS(DIST(K-1 )-DlST(K)  ) -DMI N ) 30  ,  13 , 1 3 

M  =  2 

^L  =  C 

A{  I  )=A(  I  )  +DELT(  I  )^*-CODE 

GO  TO  1  .       . 

**  TEST  FOP  CASE  3A  THROUGH  3D 

TF(ni^T(l<r-1  )-niST(X)  )21  ,26,26 

f^f^LT(  I  )=DELT(  I.)  /lO,  .   ■ 

IF(niST(K-^)-niST(i<-l  ))  73,2^,24 

A(  n=A(  I  )-DELT(  I  )*C0DE*10. 

M  =  l 

IF(DELT( I ) -DALOW) 30,30,25 


102 
o6 


10 
11 
1  2 
13 


2^ 

2^ 
24 


39 


31 


32 
?7 

30 

101 


^   :, 3  0  ♦  1 
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IF(PL-.8)51,51,52 


IF(nTST(K  l-ni.'^KK-T  )-niMIN) 

A ( T )=A( I ) +OFLT ( I )*rC0F 

IF(DIST(K-3)-niST(K-l ) )31»31 

A{ I )=Ai I )-DELT{ I )*CCDE*?. 

DFLT( I )=DFLT( I ) /lO.  .  ; 

^  =  1 

IF{DFLT(  n-DALCW)30»2  7,27 

IF(DIST !K-] )-niST(<)-DMIN) 30,30»91 

«--'  TFRMINATF  PARAMETER 

PFLT(  I  )=DFLT(  I  )^^1C. 

PRINT! 01 ,nPLT{ T ) ♦ (A ( I n ,1 1=1 ,JJ) 

FCP^;AT(2X,F15.6♦6F^?.6) 
irCNT=0 

C      «*  THIS  CHECKS  FOR  A  HIGH  LAST  VALUE  CF  THE  DISTANCE  FUNCTION 
C      ^*    IF  THIS  IS  SC  P(I)  IS  SET  EOUAL  TO  THF  PREVIOUS  VALUE 

IF(DIST{K-^)-niST(K-l))50»51,51 

DIST(<-1 )=DIST(K-3) 

A{  I  )=P1 

D(  I  )=DIST(K-1) 

1  =  1  +  1 

IF(  I-JJ-i-1  )41  »A1  ,40 

1  =  1 

M=  ] 

**  PHASE  TERMINATION  TEST 
S I  IN'  =  0  . 

DO  4?  I  I  =  ?  ,JJJ 

5UM=SUM+nA8S(n( I I-l )-D( II) ) 

IF(SUM-DTERM)43,1,1 

PRINT  110 

F0RM;.T(?X,2CHPARAMFTERS  OPTIMIZED) 

PRINT  101.CODE»{A(  I )  ,I  =  1,JJ) 

RETURN 

**  THIS  ROUTINE  WILL  INCREASE  DELT(I)  IF  IT  IS  TOO  SMALL 
Ml  =ML  +  1 

IF(ML-n  )  1  ,92,92 
D'^LT(  I  )=DFLT(  I  )*lC. 
ML  =  C 


«-« 


50 
52 

51 


AO 
41 


42 


43 

110 


91 
92 


** 


END 
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SUPPnUTINF    C(.A  .nTST»nFL.T) 

U'PLICIT    RFAL*8( A-H,C-Z) 

PIN'FN?ICN    Y(  l>n  ,n(  10  )  ,A(  10  )  ♦f^FLK  lO)  ,nTST(  "^  ) 

CC^'^'.CN    AA.  ALCW,nFLTX»D»nF,  IFRR»LKK»LLL»J  »JJ  »JJJ»KKi<;,K,I 

LKK=1 

1  ARFA=0. 

C      **  COUNTING  LOOP  FOR  TERMINATION  AFTER  C  IS  CALLED  400  TIMES     ** 
KKK=KKK+1 
816  IF{KKK-400)91»'51  ,89 

,«9  PRINT911.KKK 
911  -ORMAT(2X,?5HPROGRAM  TFRMINATED  AFTER»I4.6H  LOOPS) 
STOP 
C      **  FULERS  METHOD 

91  Z=OFLTX 
C      ^*  X        REPRESENTS  THE  INPUT  F(T) 
X=DELTX 
Y(J;=DELTX 
DO  925     I  I  =  1,JJ 
925  Y( I  I  )=w. 
C      ^*  Yd)     I-l  DERIVATIVE  OF  THE  DIFFERENTIAL  EQUATION 

2  DO  926  I  1  =  1 , JJ 

92(  Yd  I  )=Y(  I  I  )+Y{  I  I  +  l  )^!-DFLTX 

S  =  0. 

KL  =  J 

D'^  977  I  1  =  1  ♦Jj 

KL=KL-1 
927  S  =  S  +  Y( I  I ) *A( KL ) 

Y(  J)=X-S 
C      **  ERR      ABSOLUTE  VALUE  OF  ERROR  BETWEEN  INPUT  AND  OUTPUT      ** 

ERR=DABS(Y( 1 ) -X ) 
C      **  EER      MODIFIED  ERR  USED  IN  CALCULATING  AREA.   SEE  RRR  SUB   ** 

EFR=RRR(ERR»2) 
C      **  AREA     AREA  UNDER  CURVE  DEFINED  IN  RRR  FUNCTION  SUBROUTINE   ** 
76  ARFA=ARFA+EER*DELTX 

Z=Z+DELTX 

X=AINPUT(Z) 
C      **  LLL      1  PHAS=^  1  OR  ?.  2    PHASE  3.   3   PRINT  OUT  OF  RESPONSE** 

GO  TO  (210,220,400) ,LLL 

220  GO  TO  (221  ,222)  ,LKK 

221  IF(ERR-AL0W)230,230,223 

C      **  T        SETTLING  TIME  T(S) 


Al 


?^0  T=7. 
C      !^^<-  ARFA?    ST.CRES  ARFA  UP  TO  TIMF  T  IN  CA5F  T  IS  T(S)  • 

ARFA2=AREA 

L  MM  =  0  .  ' 

LKK  =  2 

GO  TO  228 
222  IF(FRR-ALCW)231 »231  »223 

231  IF(DARS{Y(2) )-.00]05)232»232»228 

232  L'^^M  =  LVM  +  1 
IF(L'^M-80'':')2»?»23  5 

23"=^  ARFA:=ARFA2 

GO  TO  241 
22:  LKK=1 

228  IF(Z-35.);  ,2»234 
234  T=2?. 
C      **  FVALUATICN  OF  OIATANCF  FUNCTION  FOR  PHASE  3 
241  DIS=DISTFN( AREA»T ) 

127  FORMAT(2X,5HAREA  » F  1  5 . 6 » 2X » 2H T  »Fl0.6»2X,  F16.9) 

PRINT127»AREA,T,A( I ) 
GO  TO  847 
2] 0  IF(2-AA  )2  »2»4 
C      *-  EVALUATION  OF  HIATANCE  FUNCTION  FOR  PHASES  1  AND  2 
4  riS=ARFA 
847  DO  867  I  I=2»4 
867  PIST( I I-l  )=DIST( IT ) 
DIST(4)=DIS 
RETURN 
C      'f^^    LISTING  OF  FINAL  RESPONSE  FOR  OPTIMAL  PARAMETERS 

400  DF=DE+DELTX 
IF(DE-1. )2»401 ,401 

401  DE=C. 

PRINT50C,Z »AREA,ERR,Y( 1  )  ; 

500  F0PVAT(6(2X,F16.9) ) 

IF(20.-Z)600,2,2 
600  RETURN 

END 

//GO.SYSIN      DD    * 
3  1   30.    .04 
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Traditionally  the  design  of  linear  control  systems  has  been  based  on 
trial  and  error  methods.   In  recent  years  considerable  interest  has  devel- 
oped in  using  computers  as  a  design  tool  to  simulate  the  control  system  and 
find  the  most  desirable  system  parameters.  This  is  done  by  defining  per- 
formance measures  in  terms  of  a  distance  function  that  relates  quality  of 
performance  to  the  least  distance  from  the  origin. 

This  paper  presents  a  FORTRAN  computer  program  that  finds  the  control 
system  parameters  that  minimize  a  distance  function  for  a  specified  topolog- 
ical configuration  and  input  f(t).   The  input  f(t)  must  reach  a  constant 
final  value  after  some  time  T^  and  never  exceed  that  value.   The  system  must 
be  described  by  an  N^*-  order  (2  N  8)  linear  differential  equation  with  con- 
stant coefficients. 

The  optimal  parameters  are  found  for  a  third  order  system  to  demonstrate 
the  program.   A  modified  step  input  is  used  with  a  performance  measure  that 
involves  settling  time  and  the  integral  of  several  forms  of  the  error  between 
the  input  and  the  output.   The  results  are  compared  and  analyzed. 

The  choice  of  the  distance  function  is  purely  subjective  and  must  be 
based  on  good  engineering  judgment.  On  this  basis  it  is  desirable  that  the 
distance  function  be  chosen  so  that  the  solution  space  is  continuous.   This 
simplifies  the  search  procedure  and  reduces  the  computation  time  required. 


